Unjamming due to local perturbations in granular packings with and without gravity 
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We investigate the unjamming response of disordered packings of frictional hard disks with the 
help of computer simulations. First, we generate jammed configurations of the disks and then force 
them to move again by local perturbations. We study the spatial distribution of the stress and 
displacement response and find long range effects of the perturbation in both cases. We record the 
penetration depth of the displacements and the critical force that is needed to make the system yield. 
These quantities are tested in two types of systems: in ideal homogeneous packings in zero gravity 
and in packings settled under gravity. The penetration depth and the critical force are sensitive to 
the interparticle friction coefficient. Qualitatively, the same nonmonotonic friction dependence is 
found both with and without gravity, however the location of the extrema are at different friction 
values. We discuss the role of the connectivity of the contact network and of the pressure gradient 
in the unjamming response. 
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I. INTRODUCTION 

Many systems including granular materials, foams and 
emulsions can flow like fluids when a high external stress 
is applied but jam into a solidHke state below a certain 
threshold of stress. In a jammed state P, [3, d, |3| the 
many-body system is trapped in a metastable configura- 
tion far from equilibrium where the constituent particles 
block each others motion. For a typical jammed granu- 
lar packing, where the thermal fluctuations are negligi- 
ble, only a sufficiently high external stress can lead to an 
unjamming transition and cause rearrangements of the 
particles. 

In this paper we study the unjamming response of 
dense disordered granular media based on computer sim- 
ulations. To achieve the unjamming transition we per- 
turb the material by generating a small local deforma- 
tion. These perturbations break the static structure of 
the packings and induce motion of particles. 

The response of granular media to local perturba- 
tions have been studied widely both in experiments 
0, S i, 0, 11 and in computer simulations % H [13, [HI- 
The majority of these studies apply small force overloads 
to study the stress response inside the bulk of material. 
In these cases, the displacements originate only from elas- 
tic deformation of the particles; the system remains in 
the jammed state. Unjamming induced by local pertur- 
bation has been also investigated experimentally. Kolb 
and co-workers 0, Q studied two-dimensional packings of 
disks under gravity and applied localized cyclic pertur- 
bations to achieve real plastic rearrangements. Based on 
the displacement field they showed that the unjamming 



response is long ranged; the magnitude of the particle dis- 
placements decays as a power law of the distance from 
the perturbation source. The exponent of the decay var- 
ied between 0.6 and 1.4 depending on the preparation of 
the system. 

Here we focus on the question of what role the in- 
terparticle friction and gravity play in the unjamming 
transition. We study the unjamming response with the 
help of the contact dynamics method jlj, [l3|, [iJl which 
handles the particles as perfectly rigid bodies. Therefore 
elastic distortion of the particles are excluded and the 
generated particle movements correspond to real plastic 
rearrangements. We analyze the response of the packing 
to local perturbations by considering the individual parti- 
cle displacements, the coarse-grained displacement field, 
and, in addition, the resistance of the system against 
the deformation. We show that the stability of the pack- 
ing against local perturbations depends strongly and in a 
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FIG. 1: Schematic picture of the two types of granular pack- 
ings confined by an external pressure bath (a) and by gravity 
(b). The dashed lines mark periodic boundaries. Some typical 
perturbations are illustrated with gray particles and arrows 
showing the direction of perturbation. 
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nontrivial way on the particle-particle friction coefficient. 
Our systems are two-dimensional disordered packings of 
disks. First, we test ideal homogeneous packings pre- 
pared by an isotropic external pressure and with fully pe- 
riodic boundary conditions [3, [l^ [Fig. [ija)] . This part 
is the full exposition and expansion of the results pre- 
sented and considered from another point of view in '15'|. 
Then we study more realistic packings that are settled 
under gravity [Fig.[TJb)]. The packings are perturbed by 
moving two adjacent particles apart in the former case 
and by shifting one particle vertically in the latter case. 

This work is organized in the following manner: Sec- 
tionllllcontains the description of the simulation method. 
In Sec. mil we present our results for the perturbation of 
homogeneous packings. The perturbation of the packings 
that are settled under gravity is investigated in Sec. IIVI 
The role of the different preparation and perturbation 
methods are discussed in Sec. [V] Finally Sec. IVII con- 
cludes the paper. 

II. SIMULATION METHOD 

We perform contact dynamics simulations on 2D gran- 
ular packings of cohesionless perfectly rigid disks. The 
numerical results of the local perturbations are presented 
for two distinct settings: (i) homogeneous random pack- 
ings confined by an external pressure in the absence of 
gravity [Fig. [Ija)] and (ii) packings confined by gravity 
[Fig.[TJb)]. In both cases, the numerical experiments con- 
sist of two steps. First we prepare static configurations of 
grains; then, we probe the packings by perturbing their 
local structure. We apply the contact dynamics method 
for both procedures. The details of the numerical 
methods are described for homogeneous and inhomoge- 
neous settings in Sees. Ill Al and III B| respectively. 

In the rest of the paper we have the following con- 
ventions. The unit of the length is set to the maximum 
grain radius. As we examine 2D systems, the disks have 
polydispersity to avoid crystallization characteristic for 
two-dimensional monodisperse ensembles. We use a uni- 
form distribution of the disk radii over the range between 
0.5 and 1. The unit of the mass is set by assuming that 
the material of the grains has unit density and the masses 
of the disks are proportional to their areas. 



A. Homogeneous random packings 

We first examine the homogeneous configurations of 
disks. Here, the acceleration of gravity is set to zero in 
order to avoid force gradients in the samples. The num- 
ber of the grains contained by the packings ranges from 
500 to 8000. As mentioned above, we first generate static 
dense random packings by compressing the initial config- 
uration of the particles into a smaller space. The com- 
paction starts with a square box filled with loose gran- 
ular gas. The disks are initially placed at random with- 



out overlaps. We apply periodic boundary conditions to 
avoid wall effects. In order to achieve homogeneous pack- 
ings, we generalized the method proposed by Andersen 
17 1 to contact dynamics (for details see flG*!). The main 
idea in this method is that, instead of using pistons, com- 
paction is achieved by imposing a constant external pres- 
sure Pcxt and let the volume of the cell evolve in time. 
In fact, the volume of the system, which acts as an addi- 
tional degree of freedom, couples with a confining pres- 
sure bath, so that the volume change is controlled by the 
difference between the external and internal pressures. 

As the size of the cell shrinks due to the difference be- 
tween Pext and the internal pressure Pin, at some point 
the grains cannot avoid touching each other anymore, 
and start building up an inner pressure to avoid inter- 
penetration. Finally all motion stops because the grains 
block further compaction. The sample is considered to 
have converged to mechanical equilibrium when further 
time evolution leads to negligible changes in the particle 
positions. At this point we have a static jammed config- 
uration under external pressure. The mechanical equilib- 
rium is achieved for each grain and the corresponding Pjn 
equals Pcxt- It is worth noting that the packing configu- 
rations depend on the friction coefficient /i. We construct 
a new packing for each value of ^ before starting the per- 
turbation process. 

The perturbation is carried out in the following way: 
We choose two adjacent particles in contact and force 
them to move apart [Fig.[IJa)]. As this case has been de- 
scribed in p^ . here only a short review is given. At 
the perturbation point we enforce the contacting sur- 
faces to open up to a small gap and determine the force 
that is needed to fulfill this constraint (critical force). 
This concept is suited very well to the contact dynamics 
method where interparticle forces are handled as con- 
straint forces, i.e., they are calculated based on con- 
straint conditions which prescribe the relative motion of 
the contact surfaces [isl . [l^ . With enforcing the opening 
of the contact, we bring the system immediately to the 
yield point where the perturbation induces sliding and/or 
opening of some contacts and initiates collective rear- 
rangements of the particles at least in the vicinity of the 
perturbation point. It is beneficial to choose small gap 
size as we are interested in the onset of motion, how the 
static structure breaks due to the perturbation. Large 
deformations, e.g., creation of new contacts, are out of 
the scope of the present study. We checked that for small 
gap sizes the displacement field (up to a constant factor) 
and the critical force become independent of the size of 
the gap. Our numerical measurements are performed in 
this gap- independent region; the size of the gap ^ is set to 
10~^ [1^. This value is far larger than the displacement 
scale 10~^^ that arises from the noise level of particle 
velocities. 

The perturbation process can be performed under two 
different boundary conditions. One can either impose the 
fixed external pressure condition in continuation of the 
assembling process, or impose the fixed volume condi- 
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tion. In the former method, the pressure of the system is 
constant and the volume is allowed to change during the 
perturbation while in the latter method, the volume is 
fixed and the pressure changes due to perturbation. This 
paper mainly contains the results of perturbation method 
with fixed pressure even though wc compare some nu- 
merical results of both methods in Sec. IIII Gl and find no 
significant differences. 



B. Packings confined by gravity 

In the inhomogeneous case, the system is settled un- 
der gravity and no pressure bath is used. Consequently 
there exists a pressure gradient in the vertical direction. 
We present simulations on model systems of = 1600 
polydisperse disks. The system is spatially periodic in the 
horizontal direction, and a one-dimensional chain of disks 
with random radii is fixed at the bottom of the box to 
provide a rough bed [see Fig. [TJb)]. The particle-particle 
and particle-base friction coefficients are the same. The 
starting configuration is a dilute granular medium which 
consists of randomly distributed nonoverlapping particles 
with zero velocities. The average initial packing fraction 
(j) ranges from 0.32 to 0.40 for different packings. Next 
the system is allowed to settle on top of the rough bed 
under the influence of the gravitational acceleration g. 
We wait until the packing relaxes into equilibrium. The 
average static packing height ranges from 35.2 ±0.1 (ap- 
proximately 23 layers of grains) for samples with small 
friction coefficient /i = 10~* to 36.9±0.1 (~ 25 layers) for 
samples with large friction coefficient /i = 10. This con- 
struction method mimics the pouring of grains through 
a sieve far away from lateral walls. 

To investigate the effect of friction and also to mea- 
sure the displacement response with better statistics, 20 
packings are constructed and perturbed for each value of 
the friction coefficient fj,. 

After that we turn to the perturbation part where we 
perturb the topmost or lowermost particles in two differ- 
ent experiments in order to test the effect of local per- 
turbation in the inhomogeneous system [Fig. [Ijb)]. In 
the first case, we choose a particle at the free surface of 
the packing and force it to move vertically downwards 
by a small displacement ^. We measure the generated 
displacements of the particle centers and also the force 
response of the system on the perturbed particle. This 
latter is the vertical component of the sum of the contact 
forces acting on the perturbed particle, which plays the 
role of the critical force. 

In the case where the system is perturbed from the 
bottom we choose a particle at the lowermost part of 
the sample and force it to move vertically upwards by a 
small upwards shift ^. The displacement pattern and the 
critical force are measured. In both cases, the magnitude 
of the displacement ^ is set to the same value as the gap 
for the homogeneous system. 



III. PERTURBATION OF HOMOGENEOUS 
PACKINGS 

In this section, we will analyze the mechanical response 
of the homogeneous packing to the local perturbation 
which we introduced in Sec. Ill Al and see how the response 
changes with the friction coefficient. The results pre- 
sented in this section belong to the system size N = 3000 
unless explicitly stated otherwise. 

After opening up a contact that is selected for the 
perturbation, we study the generated displacement field 
of particle centers and the perturbation force which is 
needed to open up the contacting surfaces at the pertur- 
bation point. This critical force -Fcrit characterizes the 
strength of the system against the local perturbation. 
Furthermore, the numerical results are compared for the 
fixed pressure and the fixed volume perturbation meth- 
ods. We close this section by reporting the results of the 
generated force and stress response fields. 

A. Displacement response field 

Our aim here is to find out how far the rearrangements 
have to penetrate into the packing as a consequence of 
the prescribed local deformation. Is there a related length 
scale? 

Our results reveal that the displacement of particle 
centers due to the single contact perturbation form a dis- 
ordered vector field. This response field can be relatively 
localized [Fig. [UJa)] or more widespread [Fig. [UJb)] de- 
pending on the location of perturbation. 

As a measure for the magnitude of the displacement 
response, we define the penetration depth d as 

N 

Y,\dM-n\ 



i=l 

where the sum runs over all particles, dj is the displace- 
ment vector of the ith particle center, n is the distance 
vector from the perturbed contact to the ith particle cen- 
ter, and n is the unit contact normal of the perturbed 
contact. 

The penetration depth 6 characterizes the size of the 
rearrangement zone in the direction of the contact nor- 
mal. S exhibits large fluctuations depending on the per- 
turbed contact. For the displacement flelds shown in 
Figs.m^a) andlDJb) the values of 5 are 7.68 and 23.31, re- 
spectively. The average penetration depth (S) , which we 
calculated based on the perturbation of 1500 randomly 
chosen contacts of the same system, is 19.8 ± 0.1. 

In order to study the average properties of the displace- 
ment fields, we perturb all contacts one by one, always 
starting with the original static packing. In each case 
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FIG. 2: (a),(b) Displacement response fields in the laboratory 
firame for two different locations of the perturbation. The re- 
sulting vector field can be relatively localized (a) or more 
widespread (b) depending on the perturbed contact, (c) Dis- 
placement response field in the contact frame averaged over 
several thousand perturbations. The system contained 3000 
disks with friction coefficient 0.5. The unit of the length, in 
which X and y are measured, is set to the maximum grain 
radius. For clarity, the magnitude of the displacements has 
been increased by a factor of 10^^ in all figures. 



the particle movements are recorded in the local contact 
frame where the perturbed contact sits in the origin and 
the X axis is chosen parallel to the contact normal, i.e., 
X indicates the direction of the separation, then we cal- 
culate the average displacement field. 

Figure[5fc) shows a smooth displacement field obtained 
by averaging over all perturbed contacts. The circular 
shape in Fig. ^c) is achieved because the original square 
shape of the system has many different orientations 
when transformed into different contact frames. Simi- 
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FIG. 3: The magnitude of the average displacements d in 
terms of the distance from the perturbed contact r. Here ^ 
stands for the gap generated at the perturbed contact. Differ- 
ent slopes correspond to different friction coefficients fi. For 
each value of friction four systems of different sizes are in- 
vestigated. The total number of particles is between 500 and 
8000. The error bars remain below 5% of the values d(r)/5 
for the whole range of r. 



lar quadrupolar structures as in Fig. were found for 
shearing amorphous systems, where localized quadrupo- 
lar rearrangements appear at the onset of plastic events 

El- 

The magnitude of the average displacement vectors d 
decays with the distance r from the perturbation point 
as shown in Fig. [DJc). In order to investigate its decay, 
we calculate d by averaging out the angle of the position. 
Figure [3] shows that the magnitude of rearrangements d 
decays as a power law of the distance r. 



d r 



(2) 



Different slopes in Fig. [3] correspond to different fric- 
tion coefficients. However, the exponent a is approxi- 
mately independent of the system size. This can be seen 
more clearly in Fig. |31[a), in which the exponent is shown 
in terms of the total number of particles N for three dif- 
ferent frictions. Our results show that a lies in the range 
0.7 — 1.4. The results of a recent experiment performed 
by moving an intruder in a system of disks Q displays 
similar power-law behavior with the same range of a. 

The power-law decay of the average displacements in- 
dicates that there is no characteristic size for the rear- 
rangement zone; instead, a decay exponent may be more 
suitable to characterize the particle movements. There- 
fore the rearrangement region is not bounded on sides by 
a penetration length and the quantity 5 may not remain 
finite for an infinitely large system. Despite of these facts, 
5 is still a useful measure of the displacements for finite 
systems. Using 5, one can still compare two displacement 
fields for the same system size. Larger 5 means a larger 
rearrangement zone. Moreover, 5 can be easily measured 
also for single perturbations where the usage of a would 
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FIG. 4: Dependence of the penetration exponent a (a) and 
the mean critical force {Fait) (b) on the system size A'^ for 
three different friction coefficients ^. 
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be troublesome. The exponent a works well for the aver- 
age displacement field but it is not a well defined quantity 
for single perturbations where the fluctuation of d{r) is 
so large that it cannot be fitted with a power-law decay. 

To investigate the role of friction, we perturb several 
packings constructed already with different friction co- 
efficients, fj,. The average penetration depth (S) has a 
strong dependence on the friction coefficient /i. ft is a 
nonmonotonic function with a sharp minimum at /i « O.I 
[Fig. O (solid circles)]. Equivalent behavior is found for 
—a as a function of /i in Ref. When the friction is 
increased starting from zero, the decrease of (6) indicates 
that the induced rearrangements become more localized. 
At /i « 0.1 the process takes a sharp turn and further 
increase of the friction leads to delocalization. 



FIG. 5: Average critical force {F^Tit) (open circles) scaled by 
the average normal contact force (-Fcont) and the average pen- 
etration depth {S) (full circles) as functions of the friction co- 
efficient /i. The vertical dashed line emphasizes that the two 
extrema are located at the same fi. 

ing that the packing constructed with friction /i = 0.1 
is the most stable packing against local perturbations. 
Moving towards the higher or lower friction coefficients, 
packings get weaker against the perturbation, the critical 
force becomes smaller, and the induced rearrangements 
become more widespread. 

To get more insight into the force response we examine 
Fcrit at every contact separately. Figure [6] reveals that 
the critical force is strongly correlated with the original 



B. Critical force 

We now turn our attention to the critical force -Fcrit- 
The results show that i^crit depends also strongly on the 
place of the perturbation. First we check its average 
properties. 

The average critical force (Fcrit) again shows strong de- 
pendence on the friction coefficient [Fig. [5] (open circles)] 
and remains approximately constant under changing the 
system size [Fig. |3Ub)]. Here, (i^crit) is scaled by the av- 
erage normal contact force (Fcont)- In (^cont), only the 
normal component of the contact forces is taken into ac- 
count and the average (• ■ •) is taken over all the contacts 
of the given packing. 

Both the critical force and the penetration depth show 
the characteristic nonmonotonic behavior (Fig. [5]) in con- 
trast to other quantities that describe the properties of 
the packing. E.g., the average coordination number 
the average contact force (i^cont), and the packing frac- 
tion exhibit smooth and monotonic functions of the fric- 
tion coefficient with plateaus for the low and high friction 
regions [l3|- -Fcrit and S behave completely differently. 
They exhibit sharp extrema at the same friction: At 
/X « 0.1 the maximum critical force and the most local- 
ized rearrangement zone is observed which has the mean- 
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FIG. 6: Each data point represents one contact in the frame 
of the critical force -Fcrit and the normal component of the 
original contact force Fcont. The four figures correspond to 
four packings constructed with different frictions. (-Fcont) is 
the average normal contact force in each packing. The dashed 
lines correspond to -Fcrit = JLont . 
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FIG. 7: Probability distribution of the critical forces Fcrit (a) 
and the critical forces normalized with respect to their mean 
in each packing i^crit/(fcrit) (b) for different friction coeffi- 
cients. The inset displays a semilogarithmic plot of proba- 
bility distribution of the normalized critical forces when av- 
eraged over all friction coefficients. The dashed line is an 
exponential fit of the tail of the curve. 



contact force. For small values of ^ a pair of contact- 
ing particles cannot resist a force of separation larger 
than the force itself that originally presses the two con- 
tact surfaces together. This can be understood well in 
the case of zero friction where the structure is isostatic 
[l9| , where the structure and the external load Pext deter- 
mine uniquely what equilibrium force is acting between 
a selected pair of contacting particles. If one pushes the 
two particles apart with a larger force, it should be com- 
pensated by a negative contact force. As the contact 
cannot exert negative forces we lose one constraint and, 
consequently, one floppy mode [l^ appears in the system 
allowing collective motion of the particles. 

Naturally, the critical force never falls below the actual 
contact force, however, for the frictional case it can get 
larger. This can be best seen for = 0.1 in Fig. [6] For 
further increase of ^ the picture becomes similar to the 
case of zero friction, i.e., the force response of a contact 
against opening is basically given by the normal compo- 
nent of the original contact force. 

It is known that contact forces in random packings 
of frictional rigid grains are not determined uniquely 
by the mechanical equilibrium and Coulomb condition 
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FIG. 8: The correlation between the critical force fcrit (full 
circles) or the normal component of the original contact force 
-fcont (open circles) and the penetration depth 5 in terms of 
the friction coefficient ^. 

[13, [l^j HH, ■ There is an ensemble of admissible force 
networks that satisfy all of these conditions in the same 
configuration of grains. The extent of the force inde- 
terminacy as a function of the friction coefficient /i was 
numerically examined in pl| for packings of rigid disks, 
where a nonmonotonic friction dependence was found 
with a maximum value at ^ « 0.1. A direct connection 
between the critical force i^ciit and the extent of the inde- 
terminacy at a given contact is established in [Tsj where 
it was found that the critical force equals the maximum 
possible contact force at the same contact taken over the 
ensemble of admissible force networks. For further details 
on the consequence of force indeterminacy the reader is 
referred to [11, HU, HI] . 

As we mentioned before, the value of the critical 
force widely changes depending on the perturbed con- 
tact. The probability distributions for Fc-dt are displayed 
in Fig. El^a) for different friction coefficients. The crit- 
ical forces are scaled in units of Fq set by the exter- 
nal pressure and by the average radius of the particles, 
Fq — 2i?avgf'cxt- Here, we use the unit Fq because, un- 
like (Fcont), it provides a fixed force scale for all sys- 
tems which is independent of the friction coefficient. Fig- 
ure El^a) shows that the probability distributions depend 
strongly on ^l. With increasing friction, probability dis- 
tributions become broader and the curves are shifted. 
The shift is nonmonotonic as expected from the behavior 
of (Fcrit). The curves move rightwards and then leftwards 
below and above fi = 0.1, respectively. 

When normalized with respect to their mean values, 
these distributions are approximately independent of the 
friction coefficient [Fig. [7jb)] and their tail can be fitted 
with an exponential decay 

P(i^crit/ (Fcrit)) OC P) /(F...) (3) 

where /3 = 2.3 ± 0.1 when averaged over all friction coef- 
ficients [see the inset of Fig. [7]^b)]. 
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FIG. 9: Comparison of the results of the fixed pressure 
(solid circles) and the fixed volume (open circles) perturba- 
tion methods, (a) The penetration depth 5 and (b) the critical 
force -Fcrit for several contacts subjected to the perturbation, 
(c) The relative change in the size of the system A (solid cir- 
cles) and the pressure P (open circles) obtained during the 
perturbation of each contact with the fixed pressure and fixed 
volume methods, respectively. Here, the friction coefficient is 
0.5. 



The tail of the curves in Fig.[71[b) for the critical forces 
is reminiscent of the tail of the probability distributions 
of the contact forces which has been studied extensively 
in the literature [H, [H, [H, [13, [H . This similarity can be 
understood well, based on the strong correlation between 
the critical force and the original contact force which is 
described in Fig. [G] 
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FIG. 10: The average magnitude of the change in the normal 
component of contact forces due to the perturbation {|AF|) 
scaled by Fo in terms of the distance from the perturbed con- 
tact r. Different curves correspond to different friction coef- 
ficients fi. 



C. Fixed pressure and fixed volume perturbations 



In Sec. Ill Al we mentioned that one can perturb the 
system by imposing either the fixed external pressure 
condition or the fixed volume condition. In this section 
we compare these two perturbation methods by applying 
them in the same test system on the same list of contacts. 
The penetration depth S and the critical force i^crit are 
shown for each perturbed contact in Figs.lHfa) and[3I^b), 
respectively. The results of the two methods are very 
similar to each other even though the imposed boundary 
conditions are basically different. 

In the fixed volume method the pressure P is allowed 
to change, while in the fixed pressure method the variable 
quantity is the volume of the system. Let us measure the 
volume change by A = AL/L, where L is the size of the 
system. Figure shows that there is a strong corre- 
lation between the variable quantities P for one method 
and A for the other. The expansion (contraction) of the 
system due to the perturbation of a contact with the 
fixed pressure method corresponds to pressure increase 
(decrease) when perturbing the same contact with the 
fixed volume method. 



Next we focus on the correlation between the critical 
force i^ciit and the penetration depth 5. Figure [5] (solid 
circles) displays that there is a weak correlation between 
these two quantities with a maximum value again around 
/i = 0.1 and it vanishes for large friction coefficients. The 
existence of the correlation means that, on average, a 
slightly larger rearrangement zone is expected for a larger 
critical force. Open circles in Fig. [5] reveal smaller cor- 
relations between the normal component of the original 
contact force -Fcont and 6. 



D. Force and stress response fields 

Finally, we investigate how other contact forces and the 
corresponding stress field is changed by the local pertur- 
bation. We measure the average magnitude of the change 
in the normal component of the contact forces (|AF|) 
which is caused by the perturbation. Figure [10] shows 
(|AF|) as a function of the distance from the perturbed 
contact. These curves, unlike the displacement response, 
do not follow a power law decay. Close to the pertur- 
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bation point, much larger (|AF|) is observed for /i = 0.1 
than for the extreme values of friction (/i = 10"^, 10). De- 
pending on the friction, {|AF|) can be changed even by 
a factor 30 in the vicinity of the perturbation point. The 
values of (|AF|) become closer to each other for different 
frictions far away from the perturbation point since the 
decay of (|AF|) becomes steeper for intermediate values 
of friction. Apparently \A.F\ goes to zero with r leading 
to tiny changes in the contact forces far away from the 
perturbation point. Interestingly, even these tiny forces 
are able to break the solid state of the packing and induce 
rearrangements of the particles (Fig. (2]). 

We calculate the average stress field caused by the per- 
turbation, measured always in the contact frame and av- 
eraged over several thousands of perturbations. We di- 
vide the system into square grid boxes of size 2. This 
corresponds to the maximum diameter or twice the min- 
imum diameter of the particles. The whole stress we 
achieved for each box is assigned to the position of the 
center of the box. The stress tensor aa/3 in each grid box 
is given [2^ by 



(a) 



1 



V 



(4) 



where V is the area of the box, is the a component 
of the force exerted on particle i by particle j, and the 
vector Hj points from the center of particle i to the center 
of j (one has to take the periodic boundary conditions 
into account and involve nearest image neighbors). The 
sum runs over all pairs of contacting particles, ^ij is a 
number between and 1 that gives what fraction of the 
line segment connecting the centers of the two particles 
is located in the box. Thus if the line segment is cut by 
grid lines the stress contribution of the contact is divided 
among the corresponding boxes. 

First we investigate the pressure which is defined as 
half of the trace of the stress tensor. Originally (before 
perturbation) the pressure is spatially constant due to the 
symmetry of the compression process. To investigate the 
effect of perturbation on the local pressure, we calculate 
the pressure change as 



AP = -tr{<T - (To) 



(5) 



where ctq is the stress tensor before the perturbation. 

Figure (TTJa) shows a logarithmic shading of the pres- 
sure change with black indicating a pressure increase, and 
white, a pressure decrease. This figure reveals that the 
pressure change decays fast with the distance from the 
perturbation point. The borders between the regions of 
positive and negative AP are also indicated in Fig. fTlT a). 
AP is positive (negative) along (perpendicular to) the 
perturbation direction. For ease of comparison Fig. lllT b) 
shows AP along and perpendicular to the x axis in the 
contact frame. 

Next we investigate the angle of shearing, which is usu- 
ally used to describe how close the material is to plastic 





FIG. 11: (a) The profile of the pressure change AP. Shading 
is logarithmic in the amplitude of the pressure change, with 
black indicating a pressure increase, and white, a pressure de- 
crease. Dots show where the sign of AP is changed, (b) AP 
along the x and y axes, (c) The maximum angle of friction 
mobilization ctmax. The x axis indicates the direction of sep- 
aration at the perturbed contact. The friction coefficient is 
0.5. 



deformation. The angle of shearing 9 is provided by the 
local stress tensor ISOll: 



arcsm I — 1 = arcsm 



P 



/<-f(a..-P)2 



P 



(6) 

When the shear stress is measured in an imaginary plane 
at some point of the material, its value depends on the 
orientation of the plane. For a given stress state, Tmax de- 
notes the maximum shear stress among all orientations. 
Dry granular media are often characterized by a critical 
angle of shearing resistance 6'crit [lOl- The material is 
expected to behave as solid until the angle of shearing 
remains below the critical angle and plastic deformation 
occurs when the critical angle is reached. 

Since the local stress is symmetric before the pertur- 
bation, the initial value of the local angle of shearing 6*0 
is approximately zero. 9 induced by the perturbation ex- 
hibits also a fast decay away from the perturbed contact. 
Figure [TlTc) shows that 9 is very small throughout the 
system. Interestingly, even the largest value 9 « 1.5° is 
far below the critical angle (the typical value of 6'crit is 
around 20° — 30°), still the perturbation is able to break 
the solid structure of the packing. 
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IV. PERTURBATION OF PACKINGS 
CONFINED BY GRAVITY 

In this section we present the numerical results of the 
local perturbations of packings confined by gravity. It it 
important to note that both preparation and perturba- 
tion steps are performed in the presence of gravity. In 
Sec. Ill Bi we explained how such static configurations are 
prepared. We also described how we perturb the topmost 
and lowermost particles in two different measurements. 
Actually, these systems are more realistic in the prepa- 
ration and perturbation methods than those investigated 
in Sec. |TlTl 

Our main aim is to study whether the nonmonotonic 
friction dependence of the mechanical response found for 
the ideal homogeneous packings is reproduced in these 
realistic systems. Here, by shifting a particle downwards 
at the free surface or upwards at the bottom of the pack- 
ing, we study the generated rearrangements of the par- 
ticle centers and also the critical force on the perturbed 
particle. 

We start our investigation with the results of the dis- 
placement response field. We find that particle move- 
ments are not bounded to a small vicinity of the pertur- 
bation point but we observe displacements even in regions 
of the system that are far from the perturbed particle. 
This indicates a long range effect similar to those found 
for homogeneous packings (Sec. IIIip and in experiments 

Here again, we characterize the size of the rearrange- 
ment zone in our finite systems with the penetration 
depth S. Similarly to Sec. IIIIl S is defined by the same 
expression [Eq. ([1])] also in the present case, only now 
denotes the distance vector from the center of the per- 
turbed particle to the center of the ith particle and n is 
pointing vertically downwards (upwards) for top pertur- 
bations (bottom perturbations). Thus 6 has the meaning 
of the vertical size of the rearrangement zone. Similar to 
the homogeneous case, S depends strongly on the per- 
turbed particle; therefore we repeat the perturbation for 
many particles to obtain the average value (<5) . 

(S) is recorded separately for top and bottom pertur- 
bations and for each value of friction. In each case, the 
average value (S) is calculated over approximately 1000 
perturbed particles. These particles are selected as fol- 
lows. We divide the width of the system into several bins 
of size roughly equal to the average diameter of the par- 
ticles. Among the particle centers located at the same 
bin, we find the highest and lowest centers and the corre- 
sponding two particles are selected for the top and bot- 
tom perturbations, respectively. We repeat this proce- 
dure for each bin. 

It has to be noted that not all the selected particles 
for bottom perturbation are taken into account in the 
calculation of (6). We exclude rattler particles that do 
not take part in the force transmission because they are 
screened by a local arch. Of course, there are no true 
rattlers in the presence of gravity, because every parti- 
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FIG. 12: The average penetration depth {S) (solid circles), 
and the average critical force (Fait) (open circles) with re- 
spect to the average effective weight (Fy) as functions of the 
friction coefficient /i, for the perturbation of topmost (a) and 
lowermost (b) particles. The vertical dashed lines are rem- 
iniscent of the position of the extrema in the homogeneous 
case. 



cle inevitably has force carrying contacts due to its own 
weight. In the case of gravity we regard a particle as a 
rattler if its upward perturbation generates no force on 
the particle, i.e. if its critical force is zero. 

To investigate the friction dependence of {S) we per- 
turb several packings constructed with different friction 
coefficients fi. Figure [T2l (solid circles) shows the average 
penetration depth as a function of fi for both top and 
bottom perturbations. It turns out that the qualitative 
behavior in both measurements is similar to the homo- 
geneous case: nonmonotonic friction dependence with a 
minimum at intermediate friction is found. However, the 
places of the minima are shifted compared to the homoge- 
neous case and also compared to each other. The possible 
explanations will be discussed in the next section. 

In fact, the bottom perturbation with gravity resem- 
bles much more the homogeneous case than the top per- 
turbation. For penetration depths in Figs. [S] and [T^b) 
the minima are quite sharp and the actual values of 6 
are almost the same concerning the minimum S and the 
values for small and large /i. For the top perturbation 
the situation is different. The penetration depth is much 
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smaller over the whole range of friction, i5 has a much 
broader minimum, and the ratio between the maximum 
and minimum values 2) is much larger than in the 
other two cases. 

For the investigation of the critical force Fcrit we first 
have to find a proper unit in which the average critical 
force is measured in order to make the results compa- 
rable. This is needed because originally the perturbed 
particles experience different average load depending on 
whether they are located in the top or in the bottom layer 
and also on magnitude of the friction. If the original load 
on a particle is larger then the critical force is expected 
to be larger as well. 

We will use the quantity (Fy) as the force unit, where 
Fy is defined for a single particle as follows: 

{c\F->a} 

Here F^ is the y (vertical) component of the contact force 
exerted on the particle at contact c, and the sum runs 
over contacts of the particle with positive Fy, i.e., over 
the supporting contacts that carry the particle against 
gravity. Similarly to {S), the average (• • •) here is taken 
over the perturbed particles and rattlers again are ex- 
cluded for bottom perturbations. 

Fy has the meaning of an effective weight of the par- 
ticle that is loaded on the supporting contacts below. It 
is constituted by the own weight of the particle plus the 
weight of other particles that is transmitted from above. 
Here we deal only with vertical components of the con- 
tact forces because the perturbation is performed in the 
vertical direction. This is in analogy with Sec. IIIII where 
the direction of the perturbation was parallel to the con- 
tact normal therefore we used the normal component of 
the original contact forces as the force unit. 

The average effective weight (Fy) is clearly different 
for the top and bottom perturbation due to the pressure 
gradient. (Fy) depends also on friction; it is an increasing 
function of /j, for both cases. For the lowermost particles, 
(Fy) ranges from 28.4 ± 0.9 for packings with small fric- 
tion coefficient fi = 10~^ to 40±1 for packings with large 
friction coefficient /i = 10. Corresponding values for the 
topmost particles are 1.25 ± 0.03 and 1.44 ± 0.04. 

We record the critical force for each perturbed parti- 
cle, where the values of Fcrit show large fluctuations. We 
determine the average (Fcrit) separately for top and bot- 
tom perturbations (rattlers are not taken into account in 
the latter case). The average critical force is displayed 
in Fig. [12] (open circles) scaled by the average effective 
weight (Fy), where the data show the dependence on fric- 
tion for both types of perturbation. We find that the non- 
monotonic behavior with quite sharp maximum, which 
was observed in Sec. IIIII is reproduced here. 

It is a common feature of Figs. [I2ja) and[T2Ub) that 
(-Fcrit) is close to the average effective weight (Fy) in the 
small friction limit and the decreasing branches of the 
curves at large frictions do not reach the same level but. 



in contrast to the homogeneous case, (Fcrit) remains con- 
siderably larger. Another difference compared to Fig. O 
is that the maxima of (Fcrit) are shifted rightwards. 

Figure [T^] reveals also some differences between the 
two types of perturbation we applied in the presence of 
gravity. The scaled critical force, e.g., is much larger in 
Fig. HHa) than in Fig. [T^ b). Interestingly, the variation 
range of the scaled critical force for the bottom pertur- 
bation is very similar to the homogeneous case. It is also 
shown in Fig.[T2]that the extrema of (Fcrit) and (S) are lo- 
cated at different values of friction for top perturbation, 
while for bottom perturbation the extrema are aligned 
similarly to Fig. [5l 

V. DISCUSSION 

In the previous two sections we applied localized per- 
turbations in a few different ways; we separated con- 
tacting particles and shifted single particles at the free 
surface or at the bottom of the system. We tested the 
critical force and the penetration depth of the perturba- 
tions both in the presence and in the absence of gravity. 
The observed behavior of these parameters was basically 
the same in all cases: They show an interesting nonmono- 
tonic dependence on the coefficient of friction (see Figs. [5] 
and [T^ . Of course, there are also some differences be- 
tween the curves presented in Figs. [H] and [T^l These dif- 
ferences may have various origins; here we discuss some 
possible causes. 



A. Connectivity 

We have used two different methods of preparation: 
"homogeneous compaction" in which the grains are com- 
pressed by a confining external pressure, and "com- 
paction by gravity" where the grains are piled due to 
gravitational acceleration. These different preparation 
methods lead to different connectivity of the packings. To 
verify this, we determine the average coordination num- 
ber z — 2Nc/N', where A'c is the total number of contacts 
and N' is the number of particles that have force carrying 
contacts. It can be seen in Fig. [T3] that the preparation 
with gravity provides a larger coordination number than 
the homogeneous compaction. The deviation is signifi- 
cant for large friction coefhcients. In the region yu > 1, 
z approximately equals the critical value 3 for homoge- 
neous compaction (open circles in Fig. I13p . This reveals 
that the structure of the packing is very close to isostatic 
[S HH HH where the contact forces are uniquely de- 
termined by the equations of mechanical equilibrium of 
the particles. The packings constructed by gravity are 
far from being isostatic in the frictional case (full circles 
in Fig. [T3|| therefore large indeterminacy of the forces is 
expected [2l| . 

Larger force indeterminacy makes the packings more 
stable against local perturbations and leads to larger crit- 
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FIG. 13: Influence of the friction on the average coordination 
number z for homogeneously confined packings (open circles) 
and packings settled under gravity (full circles). The middle 
curve (stars) corresponds to the gravity case, when the con- 
tribution of the rattler particles, that transmit no load except 
their own weight, is subtracted. 



ical forces [15|. This explains why the critical force is 



considerably larger for the right than the left side for 
both Figs. [1211 a) and [T2y b). while for the homogeneous 
case approximately the same critical forces are found in 
the small and large friction limits. This effect may also 
cause the rightward shift of the maxima of (i^crit) for the 
top and bottom perturbations in the presence of gravity. 

We note that the definition of z applied here excludes 
rattlers in zero gravity but in the presence of gravity 
all the particles are taken into account, even those that 
effectively behave as rattlers. It is important to point 
out that the difference in the connectivity we found here 
cannot be traced back entirely to the handling of rattlers. 
Even if we exclude rattlers (particles with zero critical 
force in the upward perturbation) in the calculation of z 
we cannot achieve the isostatic limit 3 for packings settled 
under gravity. This is shown in Fig. [T3] (stars) where the 
average coordination number is determined only for the 
force carrying structure, without rattlers. 

The above discussed difference in the connectivity and 
its consequences are observed only for the frictional par- 
ticles. In the zero friction limit we obtain the same coor- 
dination number z — \ for both perturbation methods. 
This is the critical coordination number for frictionless 
disks showing that inner structure of these packings is 
isostatic I3ll. 



B. Pressure gradient 

The presence of pressure gradient in the case of grav- 
ity (Sec. IIV[) makes an important difference compared to 
the homogeneous pressure (Sec. IIII|) . The perturbation 
causes displacements of the particles in a relatively large 
region. The different parts of this rearrangement zone 
can experience different local pressure. This effect is sig- 



nificant especially for top perturbations where the pres- 
sure grows proportionally with vertical distance from the 
place of the perturbation. The large relative change in 
the local pressure suppresses rearrangements in deeper 
layers which leads to considerably smaller penetration 
depths and also to larger critical forces with respect to 
the related force scale of the perturbed particles. 

One can see this overall shift of (S) and {Fa-it) for top 
perturbation in the entire region of fi [Fig. I12r a)] when 
compared to the homogeneous case (Fig. [S]). For bottom 
perturbations this effect is less important as the relative 
change in the local pressure remains small in the vicin- 
ity of the perturbed particle. This explains the smaller 
differences in {S) and (fcrit) between Figs. ISlandfT^b). 



C. Connection to force indeterminacy 



As mentioned in Sec. IIII Bl contact forces are not de- 
termined statically in packings of frictional disks. There 
exists an ensemble of force networks that solve the orig- 
inal problem, i.e., they provide mechanical equilibrium 
under the given external load and satisfy the Coulomb 
condition at every contact. We can refer to this ensemble 
as the original solution set. Due to the contact perturba- 
tion we applied in Sec. lIIII the packing finally chooses the 
force network which contains the maximum possible con- 
tact force at the perturbed contact [l^. In other words, 
the perturbation drives the system into the border of the 
original solution set. Therefore the critical force is di- 
rectly connected to the extent of indeterminacy of the 
contact forces. 

It is important to point out that the relation to the 
original force ensemble is different in the case of Sec. IIVI 
where we perturbed particles instead of contacts. Here 
the force network generated by the perturbation corre- 
sponds to a different external load and, consequently, it 
is located outside of the original solution set. Therefore 
the critical force is not related directly to the indeter- 
minacy of the contact forces in case of top and bottom 
perturbations. This seems to be a significant difference 
between contact and particle perturbations which could 
be one reason for the deviation observed for (Fcrit) be- 
tween Sees. HIT] and ITVl 



D. Other effects 

The systems constructed with homogeneous com- 
paction are close to isotropic in contrast to the pack- 
ings settled under gravity, where local stress and fabric 
are anisotropic due to the special (vertical) direction of 
the compaction. In the presence of gravity we perturbed 
the packings in this special direction which may lead to 
different response properties compared to the homoge- 
neous case. 

Moreover, for the vertical perturbations the up and 
down directions may behave differently. The direction of 
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gravity has even a structural signature in the packing. 
E.g., the number of contacts of a particle that have a 
vertical position higher than the particle center exhibits 
larger fluctuations than the number of contacts that are 
located below the center. One can mention also the par- 
ticles that are acting as rattlers when perturbed upwards 
(zero critical force) but there is no rattler behavior for 
downward perturbations where the critical force is al- 
ways positive. Therefore perturbations in the up and 
down directions may lead to a different response even if 
the perturbations are performed for the same configura- 
tion of particles inside the bulk. 

Another aspect that makes a difference between top 
and bottom perturbations is the boundary condition. 
The top of the system is not bounded, thus the displace- 
ment field generated by the perturbation is allowed to 
pass through the original boundary. This is not possi- 
ble for the bottom perturbation, where the motion of the 
surrounding particles is bounded by a rigid rough bed at 
the bottom. 



VI. CONCLUSION 

In this work we presented the numerical results of the 
measurement of mechanical response to localized pertur- 



bations. Based on contact dynamics simulations we pre- 
pared 2D static granular assemblies with and without 
gravity, then wo perturbed the systems with different 
methods to achieve the unjamming transition. Despite 
all the effects that are expected to infiuence the response 
of the packings depending on the preparation and per- 
turbation methods, the qualitative behavior seems to be 
very robust. We found that both the resistance of the 
packings against the perturbations and the penetration 
depth of the generated displacement field are sensitive to 
the interparticle friction coefficient: the surprising non- 
monotonic dependence on friction is reproduced in all 
cases that were studied in the present work. 

The nonmonotonic behavior of the critical force in the 
case of contact perturbation can be understood based on 
the nonmonotonic indeterminacy of contact forces. How- 
ever, the indeterminacy of forces seems to influence the 
critical force and the penetration depth generally. Fur- 
ther studies are needed to clarify this relationship. 
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